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ABSTRACT 

We present a study of the recent star formation of 30 Doradus in the Large Magellanic Cloud 
(LMC) using the panchromatic imaging survey Hubble Tarantula Treasury Project (HTTP). In 
this paper we focus on the stars within 20 pc of the center of the massive ionizing cluster of 30 
Doradus, NGC 2070. We recovered the star formation history by comparing deep optical and 
NIR color-magnitude diagrams (CMDs) with state-of-the-art synthetic CMDs generated with 
the latest PARSEC models, which include all stellar phases from pre-main sequence to post- 
main sequence. For the first time in this region we are able to measure the star formation using 
intermediate and low mass stars simultaneously. Our results suggest that NGC2070 experienced 
a prolonged activity. In particular, we find that the star formation in the region: i) exceeded 
the average LMG rate ~ 20 Myr ago; ii) accelerated dramatically ~ 7 Myr ago; and hi) reached 
a peak value 1-3 Myr ago. We did not find significant deviations from a Kroupa initial mass 
function down to 0.5 Mq. The average internal reddening E(B—V) is found to be between 0.3 
and 0.4 mag. 
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1. Introduction 


The Large Magellanic Gloud (LMG) harbors 
the nearest giant extragalactic Hll region ( |Ken- 
nicntt||1991 ), the Tarantula Nebula (30 Doradus). 
The central region of the nebula is dominated by 
the cluster NGG 2070, a collective of several dense 
sub-clusters ( Walborn fc Blad'^|1997 Sabbi et al. 
2012) whose light is dominated by massive OB 


stars. The most prominent of these sub-clusters 
is the very central super star cluster (SSG) R 136. 

This region has a number of characteristics that 
make it extraordinary. First of all, it displays 
an extreme rate of star formation (SF), with one 
quarter of the total massive recent (< 10 Myr) 
star formation in the LMG contained within 15' 
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from 30 Doradus (Kennicutt 1991). Moreover, 
with a stellar mass of at least 2.2 x IO^Mq con¬ 


centrated within a radius of 4.7 pc (Hunter et al 


1995), this cluster can be classified as a relatively 


low-mass clone of more distant starburst clusters, 
likely building blocks of starburst galaxies. 

Over the years, there have been debates over 
whether or not NGC 2070 is a young analog of 
old globular clusters. Recent arguments seem to 
suggest a positive answer. Using photometry [An- 
dersen et al. ( 2009[ ) found that NGC2070’s ini¬ 
tial mass function (IMF) is Salpeter-like down 
to IMq, giving the first direct evidence that 
solar-like stars are formed in a starburst cluster. 
Using multi-epoch spectroscopy of 0-type stars. 


Henault-Brunet et al. (2012) measured a velocity 
dispersion between 4 and 5 km s“^, close to the ex¬ 
pected value if it were in virial equilibrium. In this 
case the velocity dispersion would be low enough 
to allow the cluster to stay bound. 

Due to the proximity of the LMG, 30 Doradus 
can be imaged with the Hubble Space Telescope 
(HST) on scales down to ^ 0.01 pc, resolving stars 
down to the sub-solar regime. By coupling HST 
imaging with ground-based spectroscopy diagnos¬ 
tics, we know that the region has undergone con¬ 
tinuous SF activity or a superposition of multiple 
bursts of SF. Indeed, over the whole region, there 
is evidence of at least three events: i) an on-going 
off-center activity, as documented by the presence 
of embedded 0-type stars and luminous infrared 
protostars along an arc of molecular gas and warm 
dust around 30 Doradus (Rubio et al. 1992| Hy- 
land et al. 1992[ Walborn & Blades 1997[ Rubio 


et al.||19^ IWalborn et ah 1999 Brandner et al 

2001); ii) a recent event, represented by the SSG 
R 136 itself, 2 Myr old or less, as inferred from 


spectroscopy of the O stars (Massey & Hunter 


1998), 3-4 Myr old from HST imaging (Hunter 


et al. 1995); hi) an “old” event, documented by 
the presence of Hodge 301, a 25 Myr old cluster 
( Grebel fc Ghu|2QQQ ) located only 3 arcmin north¬ 
west of R 136. 


However, most of high resolution studies are 
limited to a few bands and cover only small 
patches of the 30 Doradus complex. The Hub¬ 
ble Tarantula Treasury Project (HTTP; Sabbi et 
al. 2015, in prep.), an HST survey covering a 
14' X 12' wide area around 30 Doradus at high 
resolution, fills this gap. This unique data-set. 


which combines resolution and spatial coverage, 
can be exploited in varied ways, including studies 
of massive SF, low mass SF, the IMF, the red¬ 
dening distribution, etc.. In our current paper, 
we take the opportunity to explore the SF dur¬ 
ing the last 50 Myr for the central 40 pc of 30 
Doradus, a.k.a. the starburst cluster NGG 2070, 
by comparing the observational color-magnitude 
diagrams (GMDs) with state-of-the-art synthetic 
GMDs. These simulations incorporate the lat¬ 
est (V.1.2S) set of PAdova and TRieste Stellar 


Evolution Gode (PARSEG) isochrones (see Bres- 
et al. 1(2^2 and |Tang et "ah 2014[ ), the first 


san 


theoretical library to include homogeneously all 
stellar phases from pre-main sequence (PMS) to 
post-main sequence for all masses between 0.1 
and 350 Mq. Eor the first time, we derive the 
history of the region using low and intermediate 
mass stars, either in the PMS or main sequence 
(MS) phases, simultaneously. In particular, we 
exploit magnitudes and colors of the PMS turn 
ons (hereafter TOn; see Section 1^, i.e. the GMD 
loci where the PMS phase joins to the MS. This 
“hook” has been demonstrated to be particularly 


sensitive to age (see, e.g., Stauffer 

1980| Belikov 

et al.||1998| 

Baume et al.||2003t 

Sto 

te et al.||2004t 

Mayne 2010 

Gignoni et al. 2010 

). We use both the 


optical E555W vs E555W—E775W and the near 
infrared (NIR) EllOW vs E110W-E160W GMDs 
of the HTTP sample, which offer complementary 
advantages in terms of higher spatial resolution 
(the optical) and lower reddening sensitivity (the 
NIR). 


To recover the SE rate of NGG 2070 since the 
beginning of its activity, there are three impor¬ 
tant factors that must be accounted for: differen¬ 
tial reddening, LMG field contamination, and stel¬ 
lar crowding. The first causes the GMD features 
to appear redder and more scattered, mimicking 
older populations and lowering the age resolution. 
Eield contamination mostly adds low mass stars 
to the sample, mimicking much older populations 
and a steeper IME. Einally, stellar confusion af¬ 
fects completeness and, therefore, the reachable 
look-back time. Moreover, unaccounted incom¬ 
pleteness can also mimic mass segregation. 


The paper is structured as follows. In Section 2 
we briefly describe the observations. Section 3 is 
dedicated to the physics of the TOn clock. In Sec¬ 
tion 4 we identify the general properties of the stel- 


2 



















































































lar populations in the HTTP data-set, with special 
emphasis on NGC 2070. In Section 5 we perform 
artificial star tests, the only way to take crowding 
into account, and we recover the SFH of NGC 2070 
using the synthetic CMD approach. Field contam¬ 
ination is also carefully discussed. Results are pre¬ 
sented in Section 6 and compared to the literature 
in Section 7. Section 8 compares NGC 2070 with 
other starburst clusters. Our conclusions (Section 
9) close the paper. 

Throughout the paper, for the sake of simplic¬ 
ity, we will refer to the entire HTTP data-set as 
“30 Doradus”, the starburst cluster as NGC 2070 
and its core as R 136. 


2. Observations and Photometric Reduc¬ 
tion 


We observed 30 Doradus with the HST Wide- 
Field Camera 3 (WFC3) and the HST Advanced 
Camera for Surveys (ACS) as part of proposal GO- 
12939 (PI: E. Sabbi). We built this program on an 
existing HST monochromatic survey in the F775W 
filter (GO-12499, PI: Lennon), designed to mea¬ 
sure proper motions of runaway candidates. In 
both data-sets we used the Wide Field Channel 
(WFC) of ACS in parallel with either the UVIS 
or the IR channels of WFC3 to maximize the effi¬ 
ciency of the observations. The images were taken 
with the filters F275W, F336W, F555W, F658N, 
FI low and F160W between December 2012 and 
September 2013. The survey utilized 60 orbits. 
Figure [2 shows the F775W image of 30 Doradus 
with indicated the major stellar concentrations, 
NGC 2070, Hodge 301 and NGC 2060. 


Bright and faint sources are respectively identi¬ 


fied using the packages img2xym_WFC.09x10 (An¬ 


derson & King 2006) and KS2, an evolution of 
the program described in Anderson et al. ( 2008[ ). 
A detailed description of the photometric analysis 
can be found in Sabbi et al. (2015, in preparation). 

We culled the catalog of detected objects to 
only include sources with high quality in the PSF- 
fitting, Qfit > 0.75. The final catalog contains 
^ 30, 000 detected in the filter F275W, ^ 100, 000 
in F336W, - 400,000 in F555W, - 130,000 in 
F658N, - 620,000 in F775W, 520,000 stars in 
FllOW and 570,000 in the F160W. 


3. Stellar clocks: MS Turn-Off and PMS 
Turn-On 


One of the characteristics of 30 Doradus that 
makes it a particularly interesting object is the 
high concentration of massive stars that coexist 
with a plethora of coeval intermediate and low 
mass PMS stars. This enables us to reconstruct 
the past history of 30 Doradus using the MS turn¬ 
off (MSTO^ the locus of the CMD where MS 
stars exhaust their core hydrogen, and the PMS 
TOn, where low and intermediate mass PMS stars 
ignite hydrogen in their cores. In terms of stellar 
mass, the MSTO mass is the mass of the most 
massive star still on the MS at an evolutionary 
time corresponding to the age of the cluster, while 
the TOn mass is the mass of the least massive star 
that has reached the MS at the age of the cluster. 
In the following we discuss pros and cons of the 
two clocks. 


Massive stars: The main appeal of using the 
MSTO of massive stars (M> 8Mq) stems from 
the high luminosity of such sources, which trans¬ 
lates into high quality photometry and complete 
samples. The optical and infrared spectra of these 
objects are well approximated by the Rayleigh- 
Jeans tail of a black body with temperature Te//. 
Hence, given that the spectral energy distribu¬ 
tion shape is almost unchanged as a function of 
wavelength, the optical and infrared colors are 
nearly constant (and around 0 mag). This is 
clearly visible in Figurewhere we overlaid young 
isochrones of different ages (0.5, 1, 3, 7, 15 Myr) 
on to the entire HTTP catalog (shaded grey area) 
in the UV (left panel), optical (middle panel) and 
NIR (right panel) CMDs. The adopted distance 
modulus (m—M)o and reddening E(B—V) are 18.5 
( Panagia et al.| 1991] |Schaefer||2QQ8[ Pietrzyhski 


et al. 2013) and 0.3 (chosen by eye based on the 


fit of the optical CMD), respectively. In the NIR 
CMD the upper MS (UMS) looks like a vertical 
line and isochrones of different age have MSTOs 
almost indistinguishable from one another. This 
degeneracy is attenuated in the optical and almost 
lifted in the UV CMD. Observationally, studies of 
massive stars has always been hampered by their 


^Although the term MSTO refers to the end of the MS 
phase for stars of any mass, the reference here is to stars 
more massive than 8 Mq , whose MS evolutionary times are 
shorter than 50 Myr. 
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Fig. 1.— F775W image of 30 Doradus. The stellar concentrations of NGC 2070 (continuous line cicrcle), 
Hodge 301 (dashed small circle) and NGC 2060 (dashed large circle) are also indicated. 


rarity, due to the short timescales involved in their 
evolution and the steepness of the IMF. 


On the theoretical side, models of massive stars 
are still affected by major uncertainties. In partic¬ 
ular, physical mechanisms like mass loss, rotation 
and binary evolution are not well understood (see, 
e.g., de Mink et al.||2012 ). 


Intermediate and low mass stars: For ages 
younger than 20-30 Myr, the PMS TOn is an¬ 
other valuable stellar chronometer which involves 
intermediate and low mass stars instead of massive 
star^ In analogy with the MSTO, the TOn prop¬ 
erties are directly related to the age of the stel¬ 
lar population, but with evolutionary times much 


^ Stars more massive than 6 Mq have no PMS phase at all. 


shorter than the corresponding MS times. In fact, 
the age of a cluster is equal to the time spent in 
PMS phase by its most massive star still in PMS 
phase. By definition, this star is at the TOn. 
Hence, when the intrinsic luminosity of the TOn is 
detected, it is straightforward to associate it with 
the age of the cluster. From the CMD point of 
view, the potential strength of the TOn is appar¬ 
ent from the morphology of the isochrones in Fig. 

In the optical and NIR CMDs the isochrone 
portion just before the MS has a hook and then 
is significantly flatter than the MS. The TOn is 
at the vertex of the hook, quite easy to recognize. 
For ages older than 20—30 Myr the PMS phase is 
much closer to the MS and the TOn visibility de¬ 
clines. Theory also predicts that the PMS phase 
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Fig. 2.— Stellar isochrones of the labeled ages super-imposed on the entire HTTP data-set in the UV (left 
panel), optical (middle panel) and NIR (right panel) CMDs. MSTO and TOn masses (at the age of the 
corresponding isochrones) are also indicated in the UV and optical CMD, respectively. The adopted distance 
modulus (m—M)o is 18.5 and E(B—V)=0.3. 


(recognizable as the portion of the isochrones at 
the right of the MS in Fig. itself is a valuable 
age indicator: older PMS isochrones are fainter 
and closer to the MS than younger ones. How¬ 
ever, poorly understood phenomenons like resid¬ 
ual mass accretion and magnetic fields (whose in¬ 
terplay is responsible for the appearing of irreg¬ 
ular photometric variability and UV to infrared 
excesses), and observational uncertainties like dif¬ 
ferential reddening from the circumstellar material 
can dislocate these stars (especially in the first few 
Myr) from their theoretical positions in the CMD 
(see |GouliermIs||2Q12| for a review). 


An advantage of the TOn with respect to the 
MSTO derives from the evolution at nearly con¬ 
stant luminosity of PMS stars near the MS, cor¬ 
responding to an almost horizontal track in the 
optical and NIR CMDs. This leads to a luminos- 
ity(TOn) - age relation. On the other hand, the 
MS evolution of massive stars near the MSTO is 
rather vertical (except in the UV CMD), hence age 
is not uniquely related to luminosity. Among the 
drawbacks, the intrinsic faintness of TOns com¬ 


pared to the MSTOs makes TOns prone to photo¬ 
metric errors and incompleteness, issues that are 
exacerbated in the UV CMD because older TOns 
tend to be also redder. From the theoretical side, 
the TOn visibility is intimately connected with the 
PMS evolutionary times which are still model de- 
.g., [Baraffe et al 


2009 


Hosokawa 


pendent (see, e __ 

|et al.|[2QTT| [Soderblom et al. [2014 and references 
therein). 

In this paper we aim to study the SFH of 30 
Doradus with the TOns, hence we focus our anal¬ 
ysis on the optical and NIR CMDs. As shown in 
the optical CMD (middle panel of Fig. 1^, the 
most massive star that is relevant to this study is 
around 6 Mq . 

Before closing this section, we emphasize an¬ 
other interesting feature that emerges from Fig. 
Whereas the isochrones reddened with E(B—V)= 
0.3 fit well the UMS in the optical and NIR CMDs 
(middle and right panel), the same isochrones are 
clearly too blue in the UV CMD (left panel). The 
explanation for this effect resides in the different 
populations that these CMDs trace: most of the 
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UV stars are very young 30 Doradus stars, with 
minimal contamination from field stars (the RC 
is hardly visible at F275W^ 23), while the op¬ 
tical and NIR stars can be both members and 
field stars. As a consequence, the average redden¬ 
ing of the young population in 30 Doradus must 
be higher than the adopted value E(B—V)= 0.3, 
while the average reddening of the field can be as 
low as the foreground MW reddening (E(B—V) = 
0.07; Eitzpatrick fc Savage||1984 ). 


The next Section is dedicated to a cursory CMD 
inspection of the stellar populations in the whole 
HTTP data-set. 


4. Stellar populations in the HTTP cata¬ 
log 

The CMD is the most powerful tool to recover 
the history of a resolved stellar population, since 
different parts of the CMD are populated by dif¬ 
ferent masses with different evolutionary times. In 
the following analysis we use the CMD to inves¬ 
tigate which populations are present in the whole 
HTTP catalog and the role of the NGC 2070 re¬ 
gion. Eor this task we use the optical CMD, which 
offers better spatial resolution than the NIR one. 
A more quantitative analysis will be the subject 
of Section m 

4.1. The whole HTTP sample 

The E555W vs E555W-E775W CMD of the 
entire HTTP data-set is shown in Eig. with 
overlaid isodensity contours (left panel) and stel¬ 
lar isochrones (PARSEC) of different ages (right 
panel). The first remarkable feature of this CMD 
is an extended UMS, populated by a plethora of 
intermediate and high mass stars. As the overlaid 
isochrones show, the width of the UMS is hardly 
explained by a difference in age. Moreover, given 
the average youth of these stars, plausibly much 
younger than 50 Myr, a metallicity spread is un¬ 
likely. As widely discussed in other works (see, 
e.g., ISelman et al.|p!9^ |De Marchi fc Panagia] 
2Q14| ), differential reddening is the main cause of 
this effect, although stellar rotation may also have 
a role in widening the UMS. 

To the right of the UMS, the next striking fea¬ 
ture is the very elongated red clump (RC; see also 


Section 5.1) and a broad red giant branch (RGB). 
These phases are populated by stars older than 1 
Gyr and belonging to the general field population 
of the LMG. As for the UMS width, most of the 
elongation and broadening is due to severe differ¬ 


ential reddening affecting the region (see Haschke 
et al.||2Qll[ De Marchi et al. 2Q14[ |De Marchi fc 


Panagia||2Q14). 


Einally, the lower MS shows a huge color dis¬ 
persion (more than 1 mag) at relatively bright 
magnitudes (E555W^ 23). As usual, a possible 
explanation comes from the comparison with the 
isochrones. Although differential reddening con¬ 
tributes to this effect, pushing part of the lower 
MS (LMS) to the red, the existence of short-lived 
massive stars requires that a large fraction of these 
red stars are genuine PMS stars (coeval with the 
massive stars). 

More information on the nature of these pop¬ 
ulations can be inferred from their spatial distri¬ 
bution. In fact, the spatial distribution of stars in 
different evolutionary stages yields important in¬ 
formation on the star formation processes across 
the region. Assuming a velocity dispe rsion of 21- 
27 kms“^ (measured from RGB stars; Garrera et 
al. 2011), stars older than 1 Gyr are expected to 


be diffused over scale lengths of several kpc, hence 
their distribution should be rather uniform over 
the HTTP field of view (EOV; ^200 pc). On the 
other hand, stars like those in NGG 2070, young 
and with very low velocity dispersion (4—5 km s“^; 


Henault-Brunet et al. 2012), are expected to be 


close to their birthplace. 

The right panel of Eig. [^highlights GMD re¬ 
gions of selected UMS (cyan box), RG (green box), 
low mass MS (LMS) (pink box) and PMS (blue 
box) stars. UMS and PMS stars are young objects 
a few tens of Myr old, RG stars are intermediate 
age stars (0.7 — 2 Gyr), LMS are MS stars older 
than few tens of Myr. Eigure [^ shows the corre¬ 
sponding spatial density (stars per pc^). Popu¬ 
lations grow more compact as one moves toward 
younger ages (see also Eig. 5 and 6 in Harris & 
|Zaritsky||1999| ). Other interesting features are: 

1) the UMS stars (bottom left panel in Eig. 
[^ appear very clustered. Three major concentra¬ 
tions are visible: NGC 2070, the most prominent, 
about 40 pc wide (encircled in blue), Hodge 301 
(encircled in green) and NGC 2060 (encircled 
in cyan). Lighter over-densities are also visible 
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F555W-F775W F555W-F775W 

Fig. 3.— Optical CMD of all stars in the HTTP catalog with overlaid isodensity contours (left panel) and 
PARSEC isochrones for the labeled ages (distance modulus (m—M)o and reddening E(B—V) are 18.5 and 
0.3, respectively). The dashed selections indicate sample of UMS, LMS, PMS and RC stars (see text). 


throughout the EOV. In other words, although 
most of the ongoing and recent (i.e., in the last 
50 Myr) SE is concentrated in few dense loci, a 
minor recent activity is present in the entire area; 

2) as compared to the UMS distribution, RC 
distribution (bottom-right panel) is, as expected, 
rather uniform. RC stars are a pure sample of 
field stars. Their age, older than 700 Myr, rules 
out that these objects are associated with 30 Do- 
radus, which is, at most, a few Myr old. The small 
over-densities are probably contamination of inter¬ 
mediate mass PMS stars or very reddened massive 
stars from the youngest regions. Erom a theoreti¬ 
cal point of view, the intrinsic position of the RC in 
the CMD could be affected by factors like binaries, 
differences in age and metallicity, etc. (see, e.g., 

. However, these effects are 
for the observed RC elonga¬ 
tion (up to 1.5 mag in E555W—E775W), which is 
mostly due to differential reddening. The source of 
this reddening is the gas/dust located between the 
closest and farthest RC star along the line of sight. 
More specifically, background RC stars suffer the 


Castellani et al.|[2QQQ 
insufficient to account 


highest degree of absorption, which is caused by 
the combination of Milky Way (MW) and 30 Do- 
radus extinction, while foreground stars will suffer 
only from MW (E(B—V)^ 0.07) extinction; 

3) The distributions of LMS (top-left panel) 
and PMS (top-right panel) stars appear mutually 
exclusive: the top-left corner of the EOV shows a 
paucity of PMS stars, while LMS stars are clearly 
overabundant there. The main cause for this ef¬ 
fect is reddening. Most of the LMS stars belong 
to the LMC field. This is because 30 Doradus is 
younger than few Myr, hence its MS stars are nec¬ 
essarily brighter than our LMS box of Eig. (right 
panel) (at ages <10 Myr the 30 Doradus low-mass 
stars are mainly still in PMS, as suggested by the 
isochrones in Eig. [^. Under these circumstances, 
the only way to remove LMS stars from the LMS 
box is the action of reddening (see also the map 
of Eig. [^. Doing so, these reddened MS stars 
are likely to fill the PMS box, eventually produc¬ 
ing the observed anti-correlation LMS/PMS stars 


(a similar behavior is found by Gouliermis et al. 
2006 in NGC 346; see their Eig. 4). In addi- 
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Fig. 4.— Number of stars per pc^ (see the color-bar on the right in logarithmic units) for different group of 
stars (see selection in Fig. [^: LMS (top-right panel), PMS (top-right panel), UMS (bottom-left panel), RC 
(bottom-right panel) stars. Blue, green and cyan circles highlight the regions of NGC 2070, Hodge 301 and 
NGC 2060, respectively. 


tion to this, the tendency of young PMS stars to 
be concentrated where the optical depth is higher 
exacerbates this effect. The only exceptions are 
the center of NGG 2070, where the LMS stars are 
missed because of the severe incompleteness, and 
Hodge 301, whose high concentration of LMS is 
due to its higher age (so Hodge 301’s TOn gets 
into the LMS box). 

Finally, the blue circle in Fig. is the region 
we have used to recover the SFH of NGG 2070. It 
is immediately clear from the maps that the re¬ 
gion harbors most of the UMS stars of the entire 
30 Doradus complex, as well as a remarkable con¬ 
centration of PMS stars. On the other hand the 
region is quite deficient in LMS stars, probably 
lost because of the extreme crowding conditions. 

4.2. NGC 2070 

In this Section we discuss the broad GMD 
features of the stellar population of NGG 2070. 


Our selection includes all stars within 2000 pix¬ 
els (~ 20 pc) from the cluster cente]0 Figure 
shows the corresponding GMD for the F555W 
vs F555W—F775W filters with overlaid isodensity 
contours (left panel) and PARSEG isochrones of 
1, 7 and 15 Myr (right panel). The MS con¬ 
tours of NGG 2070 differs from those of 30 Do¬ 
radus as a whole (Fig. |^. The peak density is 
around F555W~ 21, with minor peaks down to 
F555W^ 22, while the 30 Doradus GMD shows a 
smooth profile down to F555W~ 24. 

In Figure we investigate the radial distribu¬ 
tion of stars in NGC 2070. The top panels show 
CMDs of stars in concentric annuli (1-2-3) of equal 
area (the radius of the innermost circle is about 12 
pc) around the center of the cluster (from left to 
right, progressively further out from the center). 


^Throughout the paper we will refer to R 136 as the 
NGC 2070 center because it is spatially well defined. 
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Fig. 5.— CMD of NGC 2070 with isodensities overlaid (left panel) and PARSEC isochrones of 1, 7, 15 Myr 
overlaid (right panel). The adopted distance modulus (m—M)o is 18.5 and E(B—V)=0.3. 


the bottom panel shows the corresponding lumi¬ 
nosity functions (LEs). We find that: 

1) The LE of the innermost region 1 shows two 
clear peaks, located at E555W« 21 and E555W« 
24. The bright peak is well fitted by an isochrone 
of 7 Myr (see the figure), evidence of a young TOn. 
Stars brighter than this magnitude are mostly MS 
cluster members, while fainter stars are members 
only if they are on the PMS (by the definition 
of TOn, fainter members have not reached the 
MS yet). Indeed, a visual inspection of the CMD 
shows a clear color bimodality below E555W~ 21, 
which is likely due to field MS stars (mostly not- 
members) on the blue side, and PMS stars (mem¬ 
bers) on the red side. The dip after the bright peak 
is caused by the short evolutionary timescale of the 
PMS phase compared to the MS. After the dip the 
LE rises again following the IME, which increases 
at lower masses. Eventually the incompleteness 
wins, creating the second decline at E555W« 24. 
The apparent lack of lower MS stars in the re¬ 
gion 1 is probably due to the higher crowding, 
which causes more severe incompleteness than in 
the more external annuli 2 and 3. In Section [5] 


we will use the synthetic CMD approach, com¬ 
bined with artificial star experiments, to test if 
these lower MS stars are compatible with LMC 
field contamination or hide some older TOns. To 
this aim, we also need to estimate the LMC con¬ 
tamination in a reference field (see next Section). 

2) In contrast to the innermost region 1, an¬ 
nular regions 2 and 3 show a monotonic increase 
towards fainter magnitudes, with no intermediate 
peak before the final drop due to incompleteness. 
In terms of age, the lack of obvious TOns means 
that regions 2 and 3 are not dominated by as 
young stars like in region 1. Nonetheless, at the 
magnitude of the region 1 peak, region 2’s LE has 
a mild excess of stars compared to region 3, which 
would point to a poorly populated TOn. As a fur¬ 
ther proof, the region 2 CMD also shows an excess 
of PMS stars relative to region 3. 

3) Stars brighter than E ^ 18 in the region 3 
show a larger color spread than in regions 1 and 
2. The RC is clearly visible in regions 1 and 2 (at 
magnitudes 20.0—20.3 and colors 1 — 1.5), while it 
is much more dispersed in 3. All of these suggest 
a differential reddening that is higher in region 3 
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Fig. 6.— NGC 2070 radial distribution. Top panels: from left to right, stars from progressively more external 
annular regions of equal area centered onto R 136. A 7 Myr old isochrone is also shown to guide the eye. 
Lower panel: the corresponding LFs. 


than in regions 1 and 2. 


5. Recovering the SFH of NGC 2070 

The technique of recovering the SFH of spa¬ 
tially resolved populations from their CMDs (e.g. 
Tosi et aL]|1991 ) has undergone continuous refine¬ 


ment and can now provide reliable SFHs for any 
resolved population within ^ 20 Mpc (see, e.g., 
[Tolstoy et al.|2QQ^|Cignoni fc Tosi |2QlQ| and refer¬ 
ences therein). These advances have been achieved 
because of improved models of stellar evolution, 
which are now computed for fine grids of stel¬ 
lar masses and metallicities and faster multi-CPU 
computing facilities, which allow one to perform 
extensive artificial star tests on the real images 
and to fully explore a wide parameter space. 

A widely used approach consists of populat¬ 


ing a 2D array of basic synthetic CMDs gener¬ 
ated from stellar models. Each basic CMD is a 
“fuzzy” isochrone, with duration At, fixed metal- 
licity and an assumed IMF. To compare the ba¬ 
sic CMDs and the observational counterparts, the 
models are convolved with photometric errors and 
incompleteness as derived from artificial star tests 
performed on the real images. The best superpo¬ 
sition of the basic CMDs defines the SFH and the 
age-metallicity relation, and is the one that mini¬ 
mizes the residuals from observational CMD. The 
approach adopted here is described in more details 
in the Appendix. 

Despite the advances in this field, the deriva¬ 
tion of the SFH from the CMD is often affected 
by systematic errors that are difficult to assess. 
From a theoretical point of view, several stel¬ 
lar phases are still uncertain (thermally pulsing 
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asymptotic giant branch, PMS, and post-MS for 
massive stars), while observationally differential 
reddening and highly variable incompleteness can 
be hard to treat. 

Concerning NGC 2070, and the general 30 Do- 
radus region, we face three major problems in try¬ 
ing to estimate the SFH: 1) membership errors due 
to field interlopers from the LMC field that mimic 
older populations; 2) the extreme crowding con¬ 
ditions exacerbate the incompleteness, shortening 
the reachable look-back time; 3) the high level of 
differential reddening spreads and dims the CMD, 
blending together young PMS stars and older MS 
stars, which introduces further age ambiguities. 

5.1. Field contamination 

To measure the SFH of NGC 2070 we need to 
estimate the local field contamination. A typ¬ 
ical approach is to decontaminate the cluster 
by subtracting the star-counts from a reference 
field. However, given reddening and incomplete¬ 
ness variations across the whole 30 Doradus, it 
is impossible to find another direction replicat¬ 
ing the observational conditions of NGC 2070. 
A way out is to find a field that resembles as 
much as possible the LMC field we would ob¬ 
serve without NGC 2070, which means low or 
negligible SF activity in the last 50 Myr, mini¬ 
mal differential reddening and high completeness 
down to F555W^ 24. This reference LMC field 
could be then artificially corrected for the more 
severe incompleteness and photometric errors of 
NGC 2070. The resulting field would differ only 
in terms of normalization and differential red¬ 
dening from the real LMC field contaminating 
NGC 2070. Finally, reddening and normaliza¬ 
tion of this field could be tuned together with the 
SFH of NGC 2070 until an adequate match of 
NGC 2070’s CMD is found. 

As a first step we searched for low reddening 
areas around 30 Doradus. A good tracer for ex¬ 
tinction is the RC color. Figurej^shows the spatial 
distribution of RC stars, color coded according to 
the F555W—F775W color (larger and redder sym¬ 
bols correspond to redder F555W—F775W). We 
find that, at odds with the smooth distribution of 
the entire RC sample (see the bottom-right panel 
of Fig. [^, redder RC stars are very concentrated 
along filaments and arcs, most likely tracing gas 
and dust in 30 Doradus. This allows us to select 


extended regions, like the wide region just above 
NGC 2070 (see the black box in the Figure), with 
relatively blue RC stars. Because these stars can 
not all be in the foreground, these regions must 
be low-extinct ion windows in the gas/dust layers 
of 30 Doradus. Interestingly, the aforementioned 
region is also poorly populated by PMS stars (see 
top-right panel of Fig. |^, which implies minimal 
SF activity in the last 50 Myr. 

After a careful inspection of all regions with low 
extinction and PMS number, we found the best 
compromise in the region delineated by a black 
line in Fig. The corresponding CMD is shown 
in Fig. here overlaid on the 30 Doradus CMD 
(grey symbols). This sample was adopted to rep¬ 
resent our reference field. As expected, its RC 
and lower MS are much tighter than the general 
CMD, signatures of scarce differential reddening 
and young SF activity, respectively. It is worth 
noting, however, that residual differential redden¬ 
ing is still present. 

5.2. Artificial star tests 

To test the level of completeness of our pho¬ 
tometric data and to have reliable estimates of 
photometric errors, we ran extensive artificial star 
experiments. The experiments consist of adding 
“fake” sources for each of the eight pass-bands, 
modeled with the PSF used in the photometric 
analysis of the frames, onto the actual images. We 
then applied the same source detection routines 
used for our science images to the fields containing 
the combined actual images and the fake sources. 
We then determined the completeness fractions, 
defined as the ratio of recovered artificial stars to 
the number of the injected ones. We considered 
an artificial star lost if it is not recovered or if it is 
recovered being 0.75 mag brighter than its input 
magnitude. In fact, this means that it has fallen 
either on a real star brighter than the artificial 
star or one of its same brightness. Thus, we are 
not recovering or measuring the artificial star but 
a real one instead. 

To preserve the crowding conditions of the data, 
fake stars were arranged in a spatial grid such that 
the separation of the centers in each star pair was 
larger than two PSF radii. We found that a sepa¬ 
ration of 20 pixels guarantees that the probability 
of recovering a fake star is the same as it would 
be if we added only one fake star. The experiment 
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Fig. 7.— Spatial distribution of RC stars color coded according to their F555W—F775W color as given on 
the vertical bar on the right. The area delineated by a black line represents a region with low extinction 
(see text). Blue, green and cyan circles highlight the regions of NGC 2070, Hodge 301 and NGC 2060, 
respectively. 


is then repeated using a series of slightly shifted 
grids. 

As an example of our fake stars procedure, the 
top panel of Figure shows the impressive 50% 
completeness map in the F555W band for a region 
2000 pixel (~ 20 pc) away from the cluster center 
(5 million fake stars). The bottom panel shows 
the corresponding real image. The most striking 
feature is the remarkable completeness variation, 
up to 8 magnitudes, in different locations within 
the region. Regions like the broad area around 
X=10750, Y=13500 are photometrically complete 
at 50% level down to V ^ 26, whereas more 
crowded sub-regions, like the central area, are 50% 
complete only for R <20. Such variations are de¬ 
termined by the interplay of four major effects: 
photon scattering off dust in 30 Doradus, contin¬ 
uum/line nebular emission, pixel saturation and 


stellar crowding. Dust scattering and gas emis¬ 
sion raise the background flux and, in turn, make 
it more difficult to resolve stars. The net result 
is visible in the diffuse “nebulous” areas of the 
completeness map (see Figure]^. More localized 
than the dust and gas effects, pixel saturation pro¬ 
duces the visible “streaks” (charge-bleeding along 
the rows of the GGD) extending off of many bright 
stars. Finally, the extreme crowding conditions 
near R 136 are responsible for the central com¬ 
pleteness “hole” (where most of the injected fake 
stars are lost). 

5.3. Reddening and data-model compari¬ 
son 

We are not assuming any a priori reddening 
distribution. In principle we do not know where 
populations of different age are located inside 
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Fig. 8.— CMD for the reference field (black dots) overlaid on the entire 30 Doradus sample (grey dots). 


NGC 2070, so reddening is an unknown function 
of age. Naively speaking, young massive stars may 
have carved the gas with their ionizing flux, hence 
lowering the extinction. On the other hand, high 
levels of SF could be sustained only where the gas 
density is higher, which would lead to the opposite 
situation. Besides, normalization and reddening of 
the reference field are also free parameters. Fur¬ 
thermore, PMS stars may appear reddened due to 
their local circumstellar material. 

When dealing with such a large parameter 
space (age, reddening and field contamination) 
there is uncertainty on whether the best solution 
is local or global. To cope with this we combined 
a genetic algorithm with a local search procedure 
(Hybrid-Genetic Algorithm, HGA; see Appendix 
[A|for details) which is more effective to avoid local 
trapping than local search alone. In our approach 
the first step is to store a library of “basic” syn¬ 
thetic CMDs, where each GMD is a Monte Carlo 
synthetic population generated with a step-wise 
SF and reddening distribution. 

When the CMDs are generated other popula¬ 
tion parameters like metallicity, distance modu¬ 
lus (m—M)o, slope of the IMF and binary frac- 


tion q are kept fixed at 0.008 (e.g.. 

Luck et al. 

1998 

, 18.5 (e.g., Panagia et al.||1991 

r 

Schaefer 

2008 

Pietrzyhski et al.||2013 ), |Kroupa 

1 

|2001) and 


are taken from De Marchi & Panagia (2014). The 


explored ages range from now up to 50 Myr ago. 
Since older isochrones tend to be more tightly 
packed in the CMD than younger isochrones, the 
duration of each age step increases with age from 
1 Myr (at the present time) to 20 Myr (50 Myr 
ago). The reddening is allowed to vary between 0 
and 1 mag with a step of 0.05 mag. 

All “basic” synthetic CMDs (see left panel of 
Fig. 10) are then degraded using the photomet¬ 


ric errors and incompleteness as derived from the 
artificial star tests discussed in the previous Sec¬ 
tion. Once this basis is generated, any complex 
synthetic CMD can be constructed as a linear com¬ 
bination of the “basic” synthetic CMDs. In or¬ 
der to take into account field contamination, the 


CMD of the reference field (see Section 5.1) is ar¬ 
tificially reddened with steps of 0.05 mag between 


‘^Primary and secondary stars are both picked from the same 
IMF. 
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Fig. 9.— Top panel: Completeness map for NGC 2070 color-coded according to the F555W magnitude 
where the sample is 50% complete (which means that 50% of the injected artificial stars are recovered). The 
highly incomplete central region is the core of NGC 2070, R 136. Bottom panel: the real image of the region 
in the filter F555W. 
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Fig. 10.— Left panel: example of “basic” synthetic CMDs color-coded with age: blue, red, green, purple, 
violet, orange and black stars have ages in the range 0-1 Myr, 2-3 Myr, 4-5 Myr, 6-7 Myr, 8-9 Myr, 10-12 
Myr, 14-16 Myr, 18-20 Myr respectively; Right panel: example of “basic” field CMDs color-coded with 
reddening: blue, red, green stars are field stars artificially reddened with E(B—V) in the range 0-0.05 mag, 
0.35-0.40 mag and 0.70-0.75 mag respectively. 



E(B—V)=0 and E(B—V)=l to produce a com¬ 
plete basis of field CMDs. Likewise any kind of 
contamination can be simulated by linearly com¬ 
bining these “basic” field CMDs (see right panel 
of Eig. 10). The only difference is that the “ba¬ 
sic” field CMDs are corrected for the difference 
in photometric errors and incompleteness as esti¬ 
mated in the field itself and in the NGC 2070 re¬ 
gion, while the “basic” synthetic CMDs are only 
corrected for the latter. However, since in the ref¬ 
erence field completeness and photometric errors 
at E555W^ 24 are 100% and 1/20 of the error in 
sub-region Al, the first correction is almost neg¬ 
ligible. The combination of “basic” CMDs (both 
synthetic and field) which minimizes the residuals 
from the observational CMD (in terms of Poisso- 
nian likelihood) is searched with the HGA code. 
The best coefficients tell us the most likely: 1) 
star formation rate as a function of time (i.e., the 
SEH); 2) total reddening (foreground -f- internal) 


as a function of time; and 3) field contamination 
(reddening and normalization). 

Which part of the CMD is used? Not all of 
the CMD is used to recover the SEH. Our analy¬ 
sis is limited on the bright end by the saturation 
magnitude and on the faint end by the 50% com¬ 
pleteness magnitude level. Nonetheless, even with 
these conservative selections, the large variations 
of completeness with magnitude in some regions 
(see upper plot in Eig. force us to derive the 
SEH in sub-regions within which the complete¬ 
ness is reasonably uniform. In fact, our artifi¬ 
cial stars are uniformly distributed across the re¬ 
gion, whereas real stars are concentrated in clumps 
and filaments, structures covering a minor frac¬ 
tion of the total area. As a consequence, most of 
the stars will suffer a real incompleteness which 
is worse than the average measured with artifi¬ 
cial stars. The division in “iso-complete” sub- 
regions mitigates this bias. Eig. [TT] shows our se- 
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lection: sub-regions of NGC 2070 where the com¬ 
pleteness in F555W is 50% for F555W > 24 (Al), 
23 <F555W< 24 (A 2 ) and 22 <F555W< 23 (A3) 
are indicated in red, orange and blue, respectively. 
The stars found in both the F555W and F775W 
filters inside these sub-regions are used to recover 
the optical SFHs. Overall, our sub-regions trace 
the stellar density, from the lowest (Al) to the 
highest (A3). 

This analysis leaves out the very center of 
NGC 2070, a 3 pc region mostly represented by 
the SSC R 136. Here stellar crowding is so se¬ 
vere that most of the intermediate and low mass 
stars are lost, hence little, if any, information is 
available from TO ns. 

To further validate the SFH, we independently 
analized the NIR data. To do this we used all stars 
detected in both the FllOW and F160W images 
inside each sub-region (defined using the F555W 
data). In this case however the faint limiting mag¬ 
nitude was anchored to the 50% completeness level 
in FllOW. 

6. Results 
6.1. SFHs 

Figures and show the recovered SF rate 
(in M 0 yr“^pc“^; top panel) and reddening as a 
function of age (bottom panel) for sub-regions Al 
and A 2 respectively, as predicted using the op¬ 
tical (green curve) and NIR (red curve) CMDs. 
In principle our code provides the full distribu¬ 
tion of reddening as a function of age, but to 
ease the visualization we only provide the mass- 
weightec|^ average reddening as a function of age. 
The shaded bands indicate one (darker) and two 
(lighter) standard deviations. Such uncertain¬ 
ties are the quadrature sum of a statistical er¬ 
ror (obtained by bootstrapping the data and re¬ 
deriving the solutions) and systematic error (ob¬ 
tained by re-deriving the solutions with different 
age-binnings and CMD binning scheme). 

Despite the different spatial resolution and red¬ 
dening sensitivity, optical and NIR predictions are 
in good agreement. Both results predict that 
mild SF activity started throughout the whole 


NGG 2070 region 20 Myr ago, and that about 7- 
8 Myr ago the birth rate accelerated. Interestingly, 
the activity in the last 1 Myr has been relatively 
quiet compared to the average in the last 5 Myr. 

Goncerning the activity that commenced 20 
Myr ago, it is important to note that our SFHs 
have been obtained by taking into account field 
contamination, hence, unless we have been very 
unlucky in our field selectiorj^ it is reasonable to 
assume that this activity is a local extras over the 
LMG field. Nonetheless, whatever the origin, its 
significance is low (zero activity is still within 1 - 
sigma error bars). 

Focusing on the specific sub-regions, the SFH 
of region Al shows a mild bimodality, with a mi¬ 
nor peak around 5-7 Myr ago and a major peak 
1-4 Myr ago (1-3 Myr ago in the optical SFH, 2-4 
Myr ago in the NIR one). The measured redden¬ 
ing distribution anti-correlates with SF: the av¬ 
erage E(B—V) is ~ 0.6 for stars with ages older 
than 7 Myr, when the SF activity was lower, and 
^ 0.4 for stars of younger ages, when the SF was 
stronger. Optical and NIR reddening derivations 
agree within errors, both in confirming the pres¬ 
ence of notable differential reddening. The large 
oscillations of the NIR solution are mostly due to 
the low sensitivity to reddening in these filters. 
Figure shows an example of the full reddening 
solution (number of considered stars vs E(B—V)) 
for region Al in three age bins (2, 6 and 20 Myr). 
As it can be seen, the full distributions are asym¬ 
metric with a long tail at high reddening. 

The SE in sub-region A 2 (Eig. 13) is about two 
times higher than in sub-region Al. The most 
important feature in the A 2 SE is a pronounced 
peak in both optical and NIR solutions 1-3 Myr 
ago. Like in sub-region Al, the reddening of sub- 
region A 2 is anti-correlated with the SE activity. 
The average reddening values are also very similar. 

Because of the high crowding we used only the 
optical GMD to derive the SEH of the most inter¬ 
nal sub-region A3. Eig. 15 shows the result (blue) 
compared to the other two SEHs over the last 15 
Myr. The bottom panel shows the corresponding 
reddening distribution. The peak activity is al¬ 
most three times higher than in A 2 , and at odds 


Weighting with the predicted mass allows us to take into 
account that high reddening stars tend to be under-sampled 
because of the incompleteness. 


® Normalization and reddening of field stars are let to vary, 
whereas the ratio between young and old field stars de¬ 
pends on the chosen field. 
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Fig. 11.— Sub-regions of NGC 2070 where the average completeness is 50% for F555W > 24 (Al; red dots), 
23 <F555W< 24 (A2; orange dots) and 22 <F555W< 23 (A3; blue dots). 


with other sub-regions, is clearly located between 
1 and 2 Myr ago. 

The total mass assembled in the regions Al, A2, 
and A3, in the last 7 Myr, is about 2.9 x 10^ Mq, 
2.5 X 10^ Mq, and 1.1 x 10^ Mq, respectively. 

Finally, our analysis does not involve the cen¬ 
tral core of NGC 2070, where the activity is pre¬ 
sumably much higher. A rough estimate can 
be made by re-scaling the SFH of region A3 to 
match the number of stars in the magnitude range 
F555W 16-18, where the sample is fairly complete 
even in the core. In this magnitude range the 
core contains two times the stars of region A3, 
but concentrated in an area four times smaller. 
This leads to a scaled peak rate of the order 
of 2.8 X 10“^ Mq yr“^ pc“^ and a total mass of 
2.2 X 10^ Mq, in excellent agreement with Hunter 
et al. s 1995 estimate. Adding this mass estimate 


to the mass of sub-regions Al, A2 and A3, we 
get a total mass for NGC 2070 of 8.7 x 10^ M(Jj 
This value is compatible with the results of |S^ 


^ We stress that this result has been obtained with a Kroupa 


(|2QQ1|) IMF. Changing from a Kroupa to a Salpeter IMF 


below 0.5 Mq increases the total by a factor 1.6. 


man et al. (1999), who found 5.5 x 10^ Mq within 


14 pc from R 136 (using a Salpeter IMF down 
to 0.5 Mq), [Andersen et al. (2009), who found 
2.7 X 10^ Mq (using a Salpeter IMF down to 0.1 
Mq) and Bosch et al. (2001), who found 10^ Mq 
using dynamical considerations. 


6.2. Fit quality 

Figure [^compares the observed optical CMDs 
(top panels) in the three sub-regions (from left to 
right, A1-A2-A3) to synthetic CMDs (middle pan¬ 
els) generated from our best solutions. The bot¬ 
tom panels show the corresponding LFs. The red 
dashed line corresponds to the faintest limit used 
to fit the data. Figure pTj shows the same analysis 
for the NIR CMDs (regions Al and A2 only). 

Simulations for sub-regions Al and A2 show a 
good agreement with both the optical and NIR 
CMDs. Within the portion of the CMDs used to 
fit the data (see red dashed lines), observational 
and model LFs show deviations that are compat¬ 
ible to the errors (computed as square root of the 
observational count rates). This suggests that a 
Kroupa IMF, which is very similar to a Salpeter 
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Fig. 12.— Recovered SFH (top panel) and reddening distribution (bottom panel) for region Al. The optical 
and NIR solution are plotted in green and red colors, respectively. The inset panel shows the distribution of 
stars in region Al. 


IMF above 0.5 Mq, is consistent with the data. 
At fainter magnitudes our simulations systemati¬ 
cally underestimate the observational star counts 
(more in the optical than the NIR). This effect is 
minor in sub-region Al, but becomes significant in 
sub-region A2. One possible explanation for this 
mismatch is the 0.75 mag selection threshold that 
we used to reject artificial stars that are actually 
blended with real stars. Relaxing this condition 
alleviates the issue, but also increases the scatter 
in the synthetic CMDs more than what we see in 
the data. However, the impact of this effect should 
be minor in the magnitude range used for fitting. 


The simulation for sub-region A3 is in less good 
agreement with the observations. As is visible 
from the LFs (see the bottom right panel of Fig. 
16), the observed star counts systematically out¬ 


number synthetic counts even at bright magni¬ 


tudes (F555Ws< 20). The straightforward ex¬ 
planation is to attribute the discrepancy to the 
incompleteness of the region, which is too inho¬ 
mogeneous to be accurately modeled by our tests. 
However, from a physical point of view, it is also 
conceivable that the transition from PMS to MS 
for intermediate/massive stars is not well repro¬ 
duced by models. Indeed, the major discrepancy 
is in the range F555W 18-19, which corresponds to 
TOn ages younger than 1 Myr. If the PMS evolu¬ 
tionary time for these objects is overestimated (so 
more PMS stars were predicted than observed) or 
the PMS birth-line (the region in the CMD where 
stars are still embedded in their gas cocoons and, 
therefore, still invisible in optical bands) is closer 
to the MS than predicted, any attempt to simul¬ 
taneously fit PMS and MS star counts will end up 
with a bias. In this case, the deficiency of MS stars 
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Fig. 13.— Same as Fig. [T^but for region A2. 



Fig. 14.— Predicted number of stars in region A1 as a function of reddening E(B—V) for three age bins, 2 
Myr (blue), 6 Myr (green) and 20 Myr (red). 
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Fig. 15.— Optical SFHs (top panel) and reddening solutions (bottom panel) for regions A3 (blue), A2 
(purple) and A1 (pink). The inset panel shows the distribution of stars in region A3. 


could suggest that the SF rate for stars younger 
than 1 Myr might be underestimated. 

Finally, the number and colors of RC stars are 
well reproduced. This indicates that field contam¬ 
ination and reddening are correctly modeled. 


6.3. Solution robustness 


To test the robustness of the solution against 
variations in the assumed IMF, binary prescrip¬ 
tion (fraction and mass ratio) and distance, we 
re-derived the SFH of region Al using alternative 
values as an example. In Figurewe compare the 
results (shaded orange histograms) with the stan¬ 
dard solution (shaded green histogram). The top- 
left panel shows the SFH obtained using a shorter 
distance modulus, (m—M)o = 18.4, as suggested 
by some Cepheids studies (e.g. Macri et al.|2QQ6 ). 
The top-right panel shows the SFH obtained with 


a binary fraction of 60% (e.g. Sana et al. 2013) 


and mass ratio randomly drawn from a constant 
distribution between 0.5 and I. The bottom-left 
and bottom-right panels show the SFH obtained 
using an IMF exponent s, above I Mq, of 1.9 and 
2.7, respectively. 

Overall, all recovered solutions are qualitatively 
and quantitatively (within 2-sigma) consistent, 
suggesting that our findings are robust. In partic¬ 
ular, the onset of the major activity 7 Myr ago is 
unchanged. Relatively larger differences are found 
for IMF changes. The predicted rate for s= 1.9 
is significantly stronger than the 5 = 2.3 case for 
ages > 20 Myr, while for s= 2.7 it is stronger for 
ages < 5 Myr. These differences are principally 
due to the SF rate-IMF degeneracy. The lack of 
stars of lower mass (caused by the flatter IMF) is 
compensated with a higher early SF rate, whereas 
the lack of intermediate and massive stars (caused 
by the steeper IMF) is compensated with a higher 
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Fig. 16.— Observational optical CMDs (top panels) and best synthetic CMDs (middle panels) for sub- 
regions Al, A2 and A3 (from left to right). Bottom panel: data vs model LFs. The dashed line indicates 
the magnitude limit used to recover the SFH. 

recent activity. tematically underestimate and overestimate, re- 

We also found that the synthetic CMDs corre- spectively, the star-counts for F555W< 20 (the re¬ 
sponding to the s = 1.9 and s = 2.7 solutions sys- verse for F555W> 23). On the other hand, chang- 
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Fig. 17.— Same as Fig. 16 but for NIR data and sub-regions A1 and A2 only. The dashed line indicates 


the magnitude limit used to recover the SFH. 
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ing the binary fraction and mass ratio does not 
affect either the SFH or the fit quality. 

7. Comparison with Previous Studies 
7.1. SFH 


Selman et ah (1999) recovered the SFH of 


NGC 2070 using a Bayesian approach applied to 
the UBV photometry down to V = 19.2, corre¬ 
sponding to stars more massive than 20 Mq . The 
SFH was found to be dominated by three episodes, 
namely a young peak at 0 < t < 1.5 Myr, an 
intermediate-age peak at 1.5 < t < 3.5 Myr, 
and an old peak 4 < t < 6 Myr ago. These 
three bursts appear to be spatially disjoint, with 
the youngest stars concentrated towards the cen¬ 
ter while the intermediate-age stars appear spher¬ 
ically distributed over a 6 pc radius, slightly off- 
center. The observations are consistent with a SF 
that propagated inward. 


Andersen et al. (2009) obtained HST/NICMOS 


F160W band images of the central 14 pc x 14.25 
pc around R 136 and combined them with archival 


WFPC2 F555W and F814W observations (Hunter 
et al. 1995 1996). By fitting the F160W vs 


F555W—F160W CMD with stellar models above 


TMq ( [Marigo et al.||2008[) and stellar models be¬ 
low TMq ( Siess et al.||2000 ), they constrained the 
age of the low mass population to be 2 — 4 Myr 
old. 


Brandi et al. (1996) derived the age distribution 


for a region of 13 x 13 arcseconds^, at a distance 
of 4 arcseconds from R 136, using a synergy of 
NIR photometry with adaptive optics and HST 
photometry. Their investigation was limited to 
stars more massive than 12 Mq. The resulting age 
distribution (see their Fig. 13) grows from 7 Myr 
ago up to the present time with three visible bursts 
at 5-6 Myr, 3-4 Myr and 1 Myr ago. Moreover, 
they did not find red giants and red super-giants 
in their FOV, concluding that the age of R 136 
must be less than 5 Myr. 


According to Walborn & Blades (1997) the cen¬ 
tral region can be divided into a core, R 136, 2 — 3 
Myr old, a peripheral triggered population < 1 
Myr old, and a group of late-0 and early-B stars 
4-5 Myr old. 


De March! et al 


(2011) used deep HST obser¬ 
vations including Ra for a field of ^ 3 arcmin x 3 


arcmin enclosing NGC 2070. They inferred that 
a significant fraction 35%) of the PMS stars 
were formed prior to 12 Myr ago, while a similar 
fraction is younger than 4 Myr. In terms of SF, 
they found that 1 Myr ago the region was about 
30 times more active than 16 Myr ago. 

Using isochrone fitting. 


Sabbi et al. (2012) 


found that the majority of the stars in the “north¬ 
east clump”, a group of stars a few pc away from 
R 136, have ages between 2 and 5 Myr ago, while 
stars in R 136 are at most 2 Myr old. 


Except for De March! et al. (2011) and Sabbi 


et al. (2012), our data are generally much deeper 
than the others in the literature and, there¬ 
fore, more sensitive to older activity. Indeed, 
our prediction about the beginning of activity in 
NGC 2070, about 20 Myr ago, is only accessible 
with our data. Moreover, our analysis implements 
the PARSEC evolutionary models, which cover 
consistently the entire evolution from the PMS 
phase to the post-MS phase, whereas all previous 
analyses had forced to combine models for the 
PMS phase ( Siess et al.|2000 Tognelli et al.|20lf ) 
and other models for later phases. 

Despite these differences, our SEH is broadly 
consistent with the aforementioned results. Of 
particular interest is the 7 Myr epoch, when our 
solution predicts a significant SE enhancement. 
This event is generally consistent, or slightly older, 
with what is found in all other studies. Similarly 
to De March! et al. (2011) we find evidence of some 
activity prior to 10 Myr ago, whose rate is, at 
most, one order of magnitude lower than the re¬ 
cent one. 

Erom a spatial point of view, the activity shifts 
to younger ages moving from our region Al to A3, 
resembling the inward progression found by |Sel- 
man et al.||1999l This is interesting, given that 


Selman et al.||1999| and this work use complemen¬ 


tary mass ranges, above 20 Mq in the case of Sel- 
man et ^|1999 and below 6 Mq in our case. The 


inward scenario is also broadly consistent with the 


De March! & Panagia 

(2014 

), Walborn & Blades 

(|1997) and 

Sabbi et al. 

(2012 

) findings. 

De March! 

& Panagia 

(2014) found that stars younger than 


4 Myr are more concentrated towards R 136 than 
stars older than 12 Myr. These authors also found 
a remarkable lack of old activity (> 12 Myr) near 
R 136, whereas the old component is not signifi¬ 
cantly different in our sub-regions. However, the 
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Fig. 18.— Sensitivity test for the SFH in sub-region Al. Solutions for different assumptions for the IMF, 
binary population and distance (orange shaded histograms) are overlaid on the standard solution (green 
shaded histogram) . The top-left and top-right panels show the result of changing distance and binary 
fraction from 18.5 and 30% (primary and secondary masses randomly paired from the same IMF) to 18.4 
and 60% (mass ratio randomly drawn from a constant distribution between 0.5 and 1), respectively. The 
bottom-left and bottom-right show the result of changing the IMF exponent above 1 Mq from 2.3 to 1.9 and 
2.7 respectively. 


data in our sub-regions A2 and A3 are not deep 
enough to allow for a conclusive argument. 

7.2. Reddening 

Our results predict that stars younger and older 
than about 10 Myr are on average reddened by 
E(B—V)^ 0.4 ±0.05 and ^ 0.6 ±0.1, respectively. 
Because young stars are spatially more concen¬ 
trated, this also translates into a negative redden¬ 
ing gradient towards the center. This trend is in 


apparent contrast with findings of Zaritsl^ ( 1999), 
who found that the young populations (star form¬ 
ing regions few Myr old) in the LMC are more 
reddened than the older ones (> 1 Gyr). How¬ 
ever, we point out that the timescales here are 
much shorter, our “old” population is only few 
Myr older the young one. Besides, all results in¬ 
dicate a large reddening dispersion (at least 0.1 
mag). 

Looking at the literature, we find a general 
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consensus for highly variable reddening in the 
NGC 2070 neighborhood. More specifically, |De| 
March! & Panagia (2014) used UMS and RC 
stars to study extinction around NGC 2070. 
Once corrected for the foreground contribution 
(E(B—V)=0.07), our average prediction for the 
young component is very close to the peak of their 
reddening distribution obtained from UMS stars 
(see the blue histogram in their Fig. 9), while our 
average prediction for the old component is well 
within their reddening dispersion. 

Our reddening for the young population is also 
in good agreement with the results of |Selman et 


al. (1999). These authors used stars more massive 


than 20 Mq and found average extinctions in the 
range Ay = 1.1 - 1.5 (E(B-V)=0.35-0.5). 


stars (scatter that cannot be explained by obser¬ 
vational errors alone). The authors could improve 
the fit by including an age spread of 10-30 Myr in 
model clusters. 

9. Conclusions 

We have presented a detailed analysis of the 
star formation history in the starburst cluster 
NGC 2070, located in the heart of 30 Doradus in 
the EMC, using deep optical and NIR CMDs from 
the Hubble Tarantula Treasury Project. We used 
a new synthetic CMD approach combined with 
the latest Padova models (PARSEC), the first to 
cover homogeneously all stellar phases from PMS 
to post-MS. 


8. Comparison with other starburst clus¬ 
ters 

Overall, NGC 2070’s SFH is similar to that in 
NGC 346, the largest star forming region in the 
Small Magellanic Cloud (SMC), whose major ac¬ 
tivity started about 6-8 Myr ago and peaked about 
3 Myr ago before dropping to a lower level (see 
Cignoni et al. 2011). If we exclude the core of 


NGC 2070, which is ten times more active than 
NGC 346, even the peak rate of region A3 is sim¬ 
ilar in amplitude to the peak rate in NGC346 

pc“ 

Myr ago; | Cignoni et aL]|2011 ). 


(about 2x10^ Mq yr ^ pc ^ between 4 and 5 


Niederhofer et al. (2015) analysed the CMDs 


of eight young massive LMC clusters and derived 
upper limits on potential age spreads by fitting 
Gaussian profiles to their SFHs. They found age 
spreads smaller than a few Myr for the youngest 
clusters (20 - 60 Myr old), which is consistent with 
our much more detailed analysis for NGC 2070. 

Outside the Local Group, resolving older popu¬ 
lations in SSCs becomes very challenging. [Larsen 


et al. 


_ dMT] ) analysed seven young massive star 

clusters (5 - 50 Myr old) in five nearby galax¬ 
ies, located at distances of 3 - 5 Mpc. They 
found that the simulated CMDs generated with 
Padova isochrones ( Bertelli et al.|200^ ) and single 
burst SFH have problems to reproduce: 1) the ob¬ 
served separation in the CMDs between MS and 
He-burning “blue loop” stars, 2) the ratio between 
red and blue supergiants and, in some clusters, 3) 
the scatter in the luminosities of the supergiant 


In our implementation we encountered a num¬ 
ber of interesting challenges. We summarize here 
our main conclusions: 

SFR: we found that NGC 2070 experienced pro¬ 
longed activity, starting at least 7 Myr ago. We 
identify three major events in the history of this 
cluster: 

• « 20 Myr ago - This epoch demarcates the 
commencement of the first significant period 
of SF. Prior to this epoch, local activity is 
not distinguishable from the average activity 
in the LMC field; 

• 7 Myr ago - The SF accelerated throughout 
the entire region; 

• 1-3 Myr ago the activity reached a peak. In 
this time range, the SF moves from the pe¬ 
riphery to the central regions. Our inner¬ 
most region (A3) shows a maximum activity 
1-2 Myr ago; 

Stellar mass: We estimate the stellar mass of 
NGC 2070 out of 20 pc to be ^ 8.7 x 10^ Mq. 
This value is close to the median mass of Galactic 
globular clusters (8.1 x IO^Mq; 

T^). 

Reddening: Concerning the reddening distribu¬ 
tion, we find an average E(B—V)~ 0.4 mag for 
the young population (<10 Myr old), and 0.6 
mag for the old one. An explanation could be that 
only in the last few Myr the SF has been vigor¬ 
ous enough to sweep away part of the gas through 


Mandushev et al. 
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stellar winds. Another possibility is that the old 
activity took place on the far side of the NGC 2070 
nebula as seen from us, therefore those stars ex¬ 
perience most of the line of sight optical depth. 

IMF: Except the innermost few pc, where the 
incompleteness is too severe to allow firm conclu¬ 
sions, a Kroupa IMF down to 0.5 Mq is compati¬ 
ble with the data. This corroborates and extends 


the result of Andersen et al. (2009) to sub-solar 


masses. To the level we can measure low-mass 
stars can form in starburst clusters in the same 
way they form in low density environments. 
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As done previously (Cignoni et al. 2006, 2011), we parametrize the synthetic CMD as a linear super¬ 
position of basic CMDs generated with step-function-like star formation, metallicity and reddening. In 
the specific case of this paper we deal with 19 steps in age and 20 steps in reddening, for a total of 380 


parameters. To explore this wide parameter space we have combined a well tested genetic algorithm (GA), 

fl 


Pikai sl^l with a local sea rch routine. As shown in various papers (see, e.g., Ng et al.|2002 Aparicio & Hidalgo 
2009 [Small et al.|2013 ), GAs allow to find global optimum more efficiently than a local search alone. In our 
approach the synergy of the GA and a local search combines the advantages of both worlds. In the next 
sections we describe the synthetic population code and the optimization routine (GA+local search). Finally, 
the capabilities of the approach are tested with artificial data. 


A. Synthetic CMDs 


The basic synthetic GMDs(j, k) are populated through the following Monte Garlo procedure: 1) synthetic 
masses and ages are extracted from the assumed IMF and the j-th SF step, respectively; and 2) synthetic 
masses and ages are converted to absolute synthetic magnitudes and colors by using a fine grid of isochrones. 
For our calculations we used the latest (V.I.2S) PARSEC isochrones, covering the entire mass spectrum of 
O.I — 350 Mq from the PMS phase to the early-AGB phase; 3) A fraction q of synthetic stars is randomly 
chosen to have a companion. The masses of companions are extracted from the same IMF and their flux is 
added to the flux of the primaries; 4) The absolute synthetic photometry is put at the distance of the LMC 
and reddened with the k-th reddening. To produce realistic simulations, all basic CMDs are degraded with 
photometric errors and incompleteness as estimated from artificial star tests (see Section 5.2). 

To have all stellar phases well populated, the synthetic CMDs are generated with a large number of stars 
(10^). Once constructed, the basic CMDs(j, k) are binned in n bins of color and m bins of magnitude. The 
final result is a library of jxk 2D histograms and any CMD can be expressed as a linear 

combination of these 2D histograms (see Equation |A1| ). The coefficients S{j^k) that multiply each of the 
CMD^^^(jf, k) are the star formation rate at the time step j and reddening step k. The sum over j and k 
of S{j^k) X CMDm,n{j^k) provides the total star-counts predicted in the CMD bin (m,n) by the star 

formation S{j,k). 


N„.,„ = ^[S(i, k) X CMD„,„(j, k)] (Al) 

Including the reference field in the models corresponds to changing Equation |AI| into Equation |A2| 

N„,„ = ^[S(j, k) X CMD„,„(i, k)] + (A2) 

3,k 

-\-SF{k) X Fm,n{k)] 

where Sf(^) regulates the number of field stars with reddening k-th, while Fm,n{k) is the actual number of 
stars in the reference field (reddened with the k-th reddening) in the CMD bin (m,n). 


B. Best solution search 

Once the observational CMD is binned as well, the next step is to search for the combination of basic 
CMDs that minimizes the CMD residuals between data and model. Eor this task we implemented a likelihood 
distance, whose minimization is not biased by low count statistics. The combination of basic CMDs that 


® Routine developed at the High Altitude Observatory, and available in the public domain 
http: / / www.hao.ucar.edu / public/research/si / pikaia/pikaia.html 


27 














minimizes the likelihood corresponds to the most likely SFH behind the data. The uncertainty around the 
recovered best solution will be the sum in quadrature of a statistical error, obtained through a data bootstrap, 
and a systematic error, obtained by re-deriving the SFH using different age-binnings and CMD-binning. 

In our approach the likelihood is minimized with an hybrid-genetic algorithm (HGA), which combines a 
classical GA with a local search. Pure GAs are iterative probabilistic algorithms for solving a problem that 
mimic processes found in the natural biological evolution. Gompared to local search algorithms, GAs explore 
the search space in more points simultaneously, hence they are far less sensitive to the initial conditions and 
show remarkable ability to escape from local minima. Shortcomings of GAs are the weak ability for local 
exploration and the slow convergence rate. The proposed HGA aims to overcome both by alternating two 
phases: a GA, whose goal is to search for a quasi-global solution, and a local search, whose goal is to increase 
solution accuracy. The synergy of the two incorporates the exploration ability of GAs and the exploitation 
ability of local search algorithm. For this study we implemented both parallel and serial hybridization. The 
parallel one is applied to each iteration and aims at enhancing offspring likelihood by means of a local search 
before moving to the next generation. The serial one consists in applying a final local search after which the 
GA population has evolved to the region containing the global solution. 


C. Test 

To test our approach we have generated synthetic data with features resembling NGG 2070 (10000 stars 
brighter than F555W= 24). To simulate the observational conditions, this fake population is convolved with 
photometric uncertainties and incompleteness from the artificial star tests of NGG 2070. The input SFH is 
a sequence of five bursts of different ages and duration 1 Myr, while the input reddening distribution is built 
with 50% of the stars with E(B—V) between 0.20 and 0.25, and 25% between 0.25 and 0.30, 25% between 
0.15 and 0.20. Figure [T^ shows the result of the reconstruction. The input SFH and reddening distribution 
are in black; the best reconstructed counterparts in red. The recovered SFH looks like a smoothed version 
of the input SFH, due to the finite time resolution associated with the discrete nature of the data, and the 
ill-conditioned nature of mathematical inverse problems. However, overall, input functions are successfully 
recovered, with most of the differences within the uncertainties. Moreover, a bursty SFH is recovered more 
accurately if the time gap between the individual bursts is longer. 
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Fig. 19.— SFH recovery test. The top panel shows the SFH (input in black, recovered in red), bottom panel 
show the E(B—V) (input in black, recovered in red). 
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